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Summary. The application of Runge-Kutta schemes designed to enjoy a large re- 
gion of absolute stability can significantly increase the efficiency of numerical meth- 
ods for PDEs based on a method of lines approach. In this work we investigate 
the improvement in the efficiency of the time integration of relaxation schemes for 
degenerate diffusion problems, using SSP Rungc-Kutta schemes and computing the 
maximal CFL coefficients. This technique can be extended to other PDEs, linear and 
nonlinear, provided the space operator has eigenvalues with a non-zero real part. 



1 Introduction 

The integration of evolution PDE's through a method of lines approach leads 
to the solution of large systems of ODE's. Often such ODE's are stiff or 
moderately stiff; therefore the possibility of increasing the stability region of 
the time integrator can lead to a significant increase in efficiency for explicit 
or semi-implicit time-integration. 

Specifically, we consider the system of PDE's 

u t + f x (u) = Dp xx (u), (1) 

where f(u) is hyperbolic, i.e. the Jacobian of / is provided with real eigenval- 
ues and a basis of real eigenvectors for each u, while p(u) is a non decreasing 
Lipshitz continuous function, with Lipshitz constant /i. We assume the system 
has been fully discretized in space on a grid of N points Xj, j = 1, . . .N, and 
we denote with U(t) — [Ui(t), . . . , UN(t)] T the vector of the grid values of the 
numerical solution at time t. The space discretized system can be written in 
the form: 

^=L(U(t)), (2) 

leading to an autonomous system of N non linear first order ODE's. If we 
consider for instance a system of conservation laws, D = 0, the operator L 
obtained from a conservative space discretization, will be written as: 
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L(U) = --(F j+1/2 -F j _ 1/2 ), 

where is the numerical flux consistent with the physical flux f(u) in 

the usual sense of the Lax-Wendroff theorem and h is the grid spacing. 

In recent years, much research has focused on the efficient integration of 
the semidiscrete system (2). In particular, in this work we will concentrate on 
the performance of optimal Runge-Kutta schemes, characterized by a large 
stability region, introduced in [SR02]. 

Optimal Rungc-Kutta schemes are built choosing an accuracy order p and 
a number of stages s, with s > p. In principle, once s and p are fixed, the 
coefficients of the Butcher tableaux defining the Runge-Kutta schemes are 
computed maximizing in some sense the stability region and keeping as con- 
straints the fulfilment of the accuracy requirements. To compute their optimal 
schemes, Spiteri and Ruuth in [SR02] start from a strong stability condition 
which requires that the operator L be non linearly stable with respect to a 
suitable norm for a certain CFL with Forward Euler integration, namely: 

\\U n + AtL(U n )\\ < \\U n \\, VAt < At FE . 

Once this assumption is satisfied, the optimal schemes proposed in [SR02] do 
yield considerable savings in CPU time for a given accuracy. 

The idea is that an s stages Runge-Kutta scheme applied to (2) can be writ- 
ten as a convex combination of s Forward Euler steps as (see also [GST01]): 



U {1) = U n (3) 



uW=Yai k \uW+At^L(UW) (4) 

[/(«+!) = (5) 

Thus if the space discretized operator L is strongly stable for At < AtpE 
with Forward Euler time integration, than the scheme in (5) will be strongly 
stable for: 

At<XAt FE , A = min^ (6) 

Pik 

Note that (6) is only a sufficient condition for stability: clearly it may be 
possible to violate (6) and still find a stable scheme. 

The problem is that several high order space discretization operators L are 
not stable under the Forward Euler scheme, and therefore it is not possible to 
use the estimate (6) to guarantee that optimal SSP schemes will improve the 
efficiency of the time integration of (2). 

In the case of pure convection, that is for D = in (1), high order space 
discretizations of f x (u) have purely imaginary eigenvalues v up to high order 
powers of the mesh width h, that is = 0(h) p and therefore are not 
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stable under Forward Euler time integration. This is the case for instance of 
the fifth order WENO space discretization, see also [Spi], but we conjecture 
that the same holds for other widely used high order space discretizations for 
convective operators, such as ENO. 

On the other hand, when D ^ 0, the eigenvalues of the semidiscrete opera- 
tor L do have a negative real part of order 1, and therefore can be made stable 
under Forward Euler for a non zero AtFE- In this work we study the stabil- 
ity of a family of high order numerical fluxes coupled with optimal RK-SSP 
schemes for (1), in the case of pure diffusion, i.e. / = 0. In this case we find 
that several widespread space discretization schemes are stable with Forward 
Euler time integration, and therefore optimal RK-SSP schemes do yield a sig- 
nificant increase in the allowable CFL. Since in this case the eigenvalues of the 
exact operator are real, we even find that the stability estimate in (6) may be 
quite pessimistic, because it underestimates the stability region of some SSP 
schemes. In these cases, it is easy to compute numerically the maximal CFL, 
obtaining a further improvement in the efficiency of the scheme. 



2 Diffusive relaxation 



We consider high order approximations of the degenerate parabolic equation 



u t = Dp xx (u), 



(7) 



using relaxation schemes, a technique initially proposed in [JX95]. Following 
[NP00] and [CNPS06], we rewrite the system as a hyperbolic system with a 
stiff source term depending on a parameter e, which formally relaxes on the 
original parabolic equation as e — ► 0, namely: 

du dv 
dt dx 



dv o dw 
h <P Z — 

dt dx 



du 



D\ dw 

dx 



(8) 



dv _ _i 

dt dx e 



■p(u)) 



where <P 2 is a suitable positive parameter. Formally, as e — > 0, w — > p(u), 
and v — > —d x p(u) and the original equation (7) is recovered. We integrate the 
system above with an IMEX Runge-Kutta scheme [PR05]. In this fashion, the 
stiff source term is implicit and does not require restrictive stability conditions, 
while the linear convective term is explicit. 

In this work we consider only the relaxed scheme, which is obtained setting 
e = in the discretized equations. Let a^k and hi be the coefficients forming 
the Butcher tableaux of the explicit Runge-Kutta scheme in the IMEX pair. 
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The computation of the first stage value of the Runge Kutta scheme reduces 
to: 

= it" 

W W = p( u n ) . (9) 
v« = -Dd x wW 

For the following stages the first equation is 

i-l 

u {l) =u n -At^2ai, k d x v ik) . (10) 



In the other equations the convective terms are dominated by the source terms 
and thus «M and are given by 

«;W = p( u «), v« = -Dd x w^. (11) 

In this fashion, due to the particular structure of the relaxation scheme we are 
considering, only the explicit tableaux of the Implicit-Explicit Rungc-Kutta 
scheme enters the actual computation, while the coefficients of the implicit 
scheme drop out as the relaxation step is computed. For this reason we can 
apply any explicit Rungc-Kutta scheme to advance in time the solution u. 
Finally the updated solution is given by: 



u n+1 



AtJ2~hd x v^ (12) 



The overall accuracy of the scheme depends on the accuracy of the numerical 
flux used to approximate d x v^ and on the accuracy of the Runge-Kutta 
explicit scheme. 

A non linear stability analysis for the first order version of this scheme 
(upwind numerical flux on a piecewise constant space reconstruction to eval- 
uate d x v and Forward Euler in time) yields a parabolic stability restriction of 
the form: 

2h 2 1 

At<——±—, (13) 
~ fj, 1 + 2/10' v ; 

where \i is the Lipshitz continuity constant of p(u), see [CNPS06]. The 
parabolic CFL implies that to obtain a scheme of order p one should use a pth 
order accurate numerical flux for d x v and a p/2 order accurate Runge-Kutta 
method in time. 

A linear stability analysis of several higher order numerical fluxes yields 
again a parabolic CFL, of the form At < C\h 2 (l — G^tvP)//!, with constant 
C\ given by Table 1. Note that the scheme is unstable for At = C\h 2 ///, but 
for a sufficiently small h it is enough to pick At = {C\ — S)h 2 / fi for a suitable 
positive 5, i.e. the scheme is stable provided C\ is slightly decreased, see Table 
3. Note that the stability requirement becomes more strict as space accuracy 
increases, while it loosens as time accuracy increases. 
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RK1 


RK2 


RK3 


P-wise constant 


2 


2 


2.51 


P-wise linear 


0.94 


0.94 


1.18 


WEN 05 


0.79 


0.79 


1 



Table 1. CFL constant C\ for a few space reconstruction algorithms (from linear 
analysis) for standard first, second and third order RK schemes 



3 Diffusive relaxation and SSP schemes 

We start fixing a standard notation for s stages Runge-Kutta schemes. As in 
[SR02] we denote SSP(s,p) the optimal strongly stable Runge-Kutta scheme of 
orderp with s stages. We point out that SSP(1,1) is the Forward Euler scheme, 
SSP(2,2) the Heun scheme and SSP(3,3) the TVD third order Runge-Kutta 
method of [S088], which is probably the high order Runge-Kutta scheme most 
frequently used in conjunction with high order space discretizations for (1). 





s=l 


s=2 


s=3 


s=4 


s=5 


p=l 


1 


2 


3 


4 


5 


p=2 




1 


2 


3 


4 


p=3 






1 


2 


2.65 



Table 2. Improved stability coefficients for several s stages order p Runge-Kutta 
schemes (from [SR02]) 



From the first column of Table 1 we note that all the numerical fluxes 
considered are at least linearly stable with the Forward Euler scheme. Thus 
the theory of strongly stable Runge-Kutta schemes can be applied in this 
case. In particular we consider the schemes introduced in [SR02], for which 
the improved stability coefficients A of (6) can be found in Table 2. To indicate 
the fact that these coefficients are found through the theory of strongly stable 
Runge Kutta schemes and to identify to which scheme they apply, we will 
label these coefficients as X SSP (s,p). 

The coefficients A SSP (s,p) are not optimal in the case of pure diffusion. 
They can be further improved considering the actual stability region of SSP 
schemes and recalling that for the purely diffusive operator we are considering, 
only the intersection of the stability region with the real line is relevant. 

In Figure 1 we find the regions of absolute stability for some SSP schemes 
given in [SR02]. On the top left of the figure we find the stability regions of 
schemes of order 1 with 1 , 2 and 3 stages. We point out that here the improve- 
ment in the CFL constant is exactly balanced by the increased computational 
effort due to the higher number of stages: there is no gain in efficiency with 
respect to the standard Forward Euler scheme. In the top right of Figure 1 
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-6 -4 -2 



Fig. 1. Regions of absolute stability for optimal SSP schemes. On the top left: order 
1, stages 1,2,3. On the top right: order 2, stages 2,3,4. On the bottom: order 3, stages 
3,4,5. The dashed line is the Forward Euler scheme. 

there are the stability regions of schemes of order 2 with 2, 3 and 4 stages and, 
for comparison, of the Forward Euler scheme (dashed line). The SSP theory 
underestimates the effective CFL and direct inspection of the stability plot 
suggests that the stability coefficient A appearing in Table 3 can be increased 
finding the intersection of the stability curve with the real axis. Finally, on 
the bottom of Figure 1 we show the stability regions of schemes of order 3 
with 3,4 and 5 stages. Again the CFL gain of Table 2 can be improved by 
direct inspection of the graph. 

Let r] StP be the abscissa of the intersection of the stability curve with the 
negative real axis. Then the maximal gain in CFL with respect to Forward 
Euler, for a problem with real eigenvalues, is given by 

i \ \ns,p\ \vs,p\ 

Aoft[S,P — r — 

\m,i\ 2 

For several schemes in the figure it is easy to see that A OP t(s,p) > A SSP (s,p). 

The optimal CFL number for a given scheme can be found multiplying 
the coefficient G\ — S determined by the space discretization by the proper 
stability coefficient A, i.e. At < {C\ — S)X OPT (s,p)h 2 / /x. 
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stages 


CFL 


order 


N f 




stages 


CFL 


order 


N f 


2 


0.78 


4 


810 




2 


1 x 0.78 


4 


810 


3 


2 x 0.78 


4 


606 




3 


2.259x0.78 


4 


537 


4 


3 x 0.78 


4 


540 




4 


3 x 0.78 


4 


540 


SSP(3,s) + WEN05 




SSP(3,s) + WEN 05 




stages 


CFL 


order 


N f 




stages 


CFL 


order 


Nf 


3 


0.78 


4.5 


1230 




3 


1.256x0.78 


4.8 


978 


4 


2 x 0.78 


5.2 


820 




4 


2.574x0.78 


5.6 


636 


5 


2.65 x 0.78 


5.2 


770 




5 


3.106x0.78 


5.4 


660 



Table 3. Order of convergence and number of numerical flux function evaluations 
Nf (with 80 grid points) for the SSP-CFL \ SSP (s,p) (left) and the maximal CFL 
Aopt(s,p) (right) for several RK schemes. 

To measure the computational complexity of a scheme we compute the 
number Nf of numerical flux evaluations needed to reach a fixed integration 
time. Thus, for a given numerical flux, the most efficient scheme has the lowest 
value of Nf. Table 3 shows the different values of Nf obtained with CFL 
chosen according to the values of A SSP (s,p) on the left and X OPT (s,p) on the 
right. The grid spacing h is the same for all values of Nf, namely, h = 1/80. 
The table contains also data on the accuracy of the space-time scheme. The 
accuracy was evaluated using four different grids, and modifying the time step 
according to h and the chosen value of A. We don't show data for first order 
schemes since A S sp(s, 1) = A OPT (s, 1). The table shows the increased efficiency 
of the SSP(s,p) schemes with s > p (left column) and the gain obtained by a 
better estimate of A (right column). In the table, the optimal values of A are 
indicated in bold face, when they are sensitively larger than the corresponding 
X SSP (s,p). In particular, the improvement between the standard SSP(3,3) and 
the five stages SSP(5,3) with the optimal CFL is quite striking. 

The data in Table 3 refer to a linear diffusion problem, although we find 
analogous results for a nonlinear degenerate diffusion equation. We performed 
tests on the self-similar Barenblatt solution of the porous media equation. In 
this case the order of accuracy is limited by the non regularity of the solution, 
but we find that the errors with respect to the exact solution slightly decrease 
using the SSP(s,p) with s > p and the optimal A. 

4 Final remarks 

We have shown that the theory of [SR02] can be applied to diffusion equations 
in the relaxation framework to improve the efficiency of the time integration. 
Moreover the fact that the eigenvalues of the semidiscrete operator are real 
numbers permits to further improve the efficiency of the schemes. We ex- 
pect that analogous results can be obtained for convection-diffusion operators, 
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thanks to the nonzero real part of the eigenvalues of the discrete operators. 
This allows to achieve stability under Forward Euler integration. 

We also note that the same framework can be applied to other space 
discretizations besides the numerical fluxes obtained via relaxation schemes: 
the key factor is the localization of the eigenvalues of the differential operator, 
which must have a non zero real part. On the other hand, the plots of absolute 
stability regions show that Runge-Kutta schemes with number of stages s > p 
can be built with improved stability regions, notwithstanding stability under 
Forward Euler. These schemes can be applied to semidiscrete operators even 
in the convective regime, but their stability conditions cannot be derived from 
the behaviour of the operator under the Forward Euler scheme. 
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